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These lectures are for novices to lattice QCD. They introduce a set of simple ideas 
and numerical techniques that can be implemented in a short period of time and 
that are capabable of generating nontrivial, nonperturbative results from lattice 
QCD. The simplest of these calculations can be completed on a standard worksta- 
tion or high-end personal computer. 



1 Introduction 

These lectures are for novices who are interested in learning how to do lattice 
QCD simulations. My intent is to describe in detail everything that one needs 
to know in order to create and run a simple lattice QCD simulation. My 
focus here is not on lengthy derivations or detailed comparisons of algorithms. 
Rather I want to introduce a set of simple ideas and techniques that one 
can implement in a relatively short time and that are capable of generating 
nontrivial results from lattice QCD. 

We begin in Section 2 with simple one-dimensional quantum mechanics. 
Most of the simulation techniques can be applied to ordinary quantum me- 
chanics, and the simulations require only seconds or minutes of computer time, 
rather than the hours or days needed for QCD simulations. Consequently such 
applications are ideal for learning the simulation technology. We broaden the 
discussion to include quantum field theories in Section 3. Here we discuss the 
theoretical techniques needed to obtain accurate results on the relatively coarse 
lattices well suited to smaller computers. Finally, in Section 4, we adapt our 
techniques for simulations of gluon dynamics in QCD. 

2 Numerical Path Integrals 

2.1 Discretizing the Path Integral 

We begin with numerical techniques for evaluating path integrals. Recall what 
path integrals tell us. In one-dimensional quantum mechanics, for example, the 
evolution of a position eigenstate |x;) from time t{ to time tf can be computed 
using a path integral 1 : 

-H(«f-ti) u.\ - / tw+'U-SM 



x { \e-"^- tl > \aa) = / Vx(t)e'^ xi . (1) 
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Here the / Vx{t) designates a sum over all possible particle paths 

{x(t) for t = U -» ti} (2) 

with 

x(U) — Xi x(t{) — X{. (3) 

The hamiltonian is H, and S[x] is the classical action, 

S[x] = / dtL(x,x) = dt 
Jti Jt, 

evaluated for each path x(t). There arc no i's in these formulas because we 
are using "euclidean" path integrals. These are derived in the same way as 
standard path integrals but with it — ► t. Euclidean path integrals are much 
better for numerical work since the integrands do not oscillate wildly in sign. 

Knowledge of the propagator, Eq. (1), as a function of Xi, ti, Xf , tf gives us 
complete information about the quantum theory. For example, we can easily 
determine the groundstate energy and wavefunction. Setting 

X; = X{ = X tf — t[ = T, (5) 

the propagator can be rewritten 

(x\ e~" T \x) = J2 {A En) c- E " T (E n | x) (6) 

n 

where \E n ) is the energy eigenstate with eigenvalue E n . The sum is dominated 
by the lowest-energy states when T is large, because of the exponentials, and 
in the limit of very large T only the groundstate, |-Eb), contributes: 

(x\e-" T \x) T -=5° e-^|(,|£ >| 2 . (7) 
We extract the groundstate energy E by integrating over x, 

J dx (x\e-" T \x) T -=T e- BoT , (8) 

and then, going back to the previous equation, we determine the groundstate 
wavefunction ipE ( x ) = ( x \ Eq). 

Our goal therefore is to develop a numerical procedure for evaluating the 
propagator using a path integral. There are two issues we must address. First 
we must find a way to represent an arbitrary particle path {x(t), U < t < tf} 



mx(ty 



+ V(x(t)) 



(4) 



2 



in the computer. A path is specified by a function x(t) which, in principle, 
can be infinitely complex (and therefore too much for any computer). We 
approximate this function by specifying x(t) only at the nodes or sites on a 
discretized t axis: 

tj =ti+ j a for j = 0, 1 . . . N (9) 
where a is the grid spacing, 

t{ ~ *i nn\ 

ttE — ( 10 ) 

Then a path is described by a vector of numbers, 

x = {x(t ),x{t 1 )...x{t N )}. (11) 

It is common practice to refer to such a path as a "configuration" . The integral 
over all paths in this approximation becomes an ordinary integral over all 
possible values for each of the a;(tj)'s: that is, 

/ Vx{t) -> A [ dx 1 dx 2 ...dx N - 1 (12) 

J J — oo 

where we have adopted the notation Xj = x(tj). We don't integrate over the 
endpoints since they are held fixed; for example, for boundary conditions (5), 

xo = x n = x. (13) 

We won't need the normalization factor A for most of our work, but for our 
one-dimensional problem it is 1 

The second issue we must address concerns the evaluation of the action 
given only a discretized path {xj}. Focusing just on the contribution from 
tj <t<tj+\, the obvious approximation is 



rtj+i 

/ dtL 

Jtj 



m ( x j+1 - Xj\ 2 1 



+ ^(V(x j+1 ) + V( Xj )) 



(15) 



With this approximation, our numerical representation of the path integral is 
complete, and we have an approximate expression for the quantum mechanical 
propagator: for example, 

/oo 
dx 1 ...dx N - 1 e- s, *W (16) 
-oo 



3 



where 

JV-l 

Si at [x] = - x i? + aV ( x j)} > ( 17 ) 

3=0 

xo = xn = x, and a — T/N. We have reduced quantum mechanics to a 
problem in numerical integration. 

One might worry about approximating x with — Xj)/a in our formula 
for the lattice action Sut [x] . It is not obvious that this is a good approximation 
given that Xj + \ — Xj can be arbitrarily large in our path integral; that is, paths 
can be arbitrarily rough. While not so important for our one-dimensional 
problem, this becomes a crucial issue for four-dimensional field theories. It is 
dealt with using renormalization theory, which we discuss in later sections. 

Exercise: Set a = 1/2 and N = 8 in approximate formula (16) for the propagator, 
and integrate the right-hand side numerically. The seven dimensional integral 
that results is easily evaluated using standard routines, such as vegas 2 . Do 
this first for the one- dimensional harmonic-oscillator potential 

x 2 

V(x) = — with m = 1. (18) 

Evaluate the propagator for several values of xo = a; at = x, and compare your 
results with those of standard quantum mechanics: 

( a; |e-" T | a; )«|( S |£;o>| 2 e- BoT (19) 

where Eo = 1/2 and 

<*l E o) = e -^T- (20) 
Extract the energy and wavefunction from your numerical result. Repeat this 
exercise for V(x) — x 4 /2. (If you wish, you may restrict x integrations to the 
region — 5 — > 5 rather than — oo — > oo; this has negligible effect on the results 
of this exercise) . 

My results for the harmonic oscillator case are shown in Fig. 1. 



2.2 Monte Carlo Evaluation of Path Integrals 

Our analysis in the previous section focused on the groundstate. In quantum 
field theory, where the groundstate is the vacuum, we are generally interested in 
excited states. To analyze excited states using path integrals, we interrupt the 
propagation of the groundstate by introducing new operators at intermediate 
times. Consider, for example, the quantity 
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{x\e~ HT \x) 



0.05 




Figure 1: Euclidean harmonic-oscillator propagator at large time (T = 4) computed exactly 
(dotted line), and computed using numerical integration to evaluate the discretized path 
integral (points) . The path integral was approximated by an 8 dimensional integral which was 
evaluated numerically, using vegas, at the points indicated. The exact result is approximated 
by the square of the ground state wavefunction multiplied by exp(— EqT). 



where now we integrate over all x\ = Xf — x as well as the intermediate x(t)'s. 
This quantity is a weighted average of x(t 2 )x(ti) over all paths, with weight 
exp(— S[x]). The numerator on the right-hand side equals, in quantum me- 
chanics, 

J dx (x\ e-^-hlie-^-^ie-^- 41 ' \x) . (22) 
Setting T = tf — U and t = t% — t\ we can rewrite the full expression as 

{{x(kHh)}) . E«-"<^-*»*l*.) . (23) 

If T » t and large, then the groundstate |£o) dominates the sums and 

G(t) = ((xfaMh))) -» (E \xc-^- E ^x\E ). (24) 

In our harmonic oscillator example, the state propagating between the two x's 
cannot be \E ) since x switches the parity of the state. Thus if we now make 
t large (but still < T) 

G(t) tl ^ c KSo|i|Si)| 2 e- (Bl - So)t (25) 
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where \Ei) is the first excited state. Consequently we can extract the first 
excitation energy from the large-i dependence of G(t), 

\og(G(t)/G(t + a)) - (Eh. - Eo)a, (26) 

and then, going back to G(t), we can determine the quantum mechanical tran- 
sition matrix element (Eq\ x\E{). 

In principle, path integral averages ((r[x])) of arbitrary functional Y\x] 
can be used to compute any physical property of the excited states in the 
quantum theory. Also we note in passing that 



Ee- £ " T (E n \T[x]\E n ) 



«TN» = ^ ^ r ' ( 27 ) 



becomes a (quantum mechanical) thermal average if we replace 

T - 13 = l/k B T tcmp . (28) 

Thus any computer code designed to compute path integral averages can be 
used for thermal physics as well. Here we focus on the zero-temperature limit 
of large T. 

We could evaluate the path integrals in ((r[x])) using a standard multi- 
dimensional integration code like vegas, at least for one-dimensional systems. 
Here, instead, we employ a more generally useful Monte Carlo procedure. Not- 
ing that 

r 1XX [Vx(t)r[x]e- s W 
{ ™ = fVx\t)e-sW ■ ™ 

is a weighted average over paths with weight exp(— S[x}), we generate a large 
number, N c {, of random paths or configurations, 

x^ = {xi a) x[ a) ...x^l,} a = l,2...JV rf , (30) 

on our grid in such a way that the probability P[a;( Q )] for obtaining any par- 
ticular path x< Q > is 

P[x^\ <xe- s ^ xM l (31) 

Then an unweighted average of T[x] over this set of paths approximates the 
weighted average over uniformly distributed paths: 

1 Ncl 

((T[x]))«T=— ^r[^)]. (32) 

a—1 



r is our "Monte Carlo estimator" for ((r[x])) on our lattice. Of course the 
estimate will never be exact since the number of paths N c { will never be infinite. 
The Monte Carlo uncertainty tJp in our estimate is a potential source of error; 
it is estimated in the usual fashion 3 : 



(33) 




This becomes 



, _ «r 2 » - «r»' 



for large N c f. Since the numerator in this expression is independent of N c t 
(in principle, it can be determined directly from quantum mechanics), the 
statistical uncertainties vanish as l/y/N c { when iV c f increases. 

We need some sort of specialized random- vector generator to create our set 
of random paths x^ with probability (31). Possibly the simplest procedure, 
though not always the best, is the Metropolis Algorithm 4 . In this procedure, 
we start with an arbitrary path x 1 - ) and modify it by visiting each of the 
sites on the lattice, and randomizing the Xj's at those sites, one at a time, in 
a particular fashion that is described below. In this way we generate a new 
random path from the old one: x^ — ► x^ . This is called "updating" the 
path. Applying the algorithm to a^ 1 ) we generate path x^ 2 \ and so on until we 
have iV c f random paths. This set of random paths has the correct distribution 
if N c f is sufficiently large. 

The algorithm for randomizing Xj at the j th site is: 

• generate a random number (, with probability uniformly distributed be- 
tween — e and e for some constant e; 

• replace Xj — > Xj + £ and compute the change AS in the action caused by 
this replacement (generally only a few terms in the lattice action involve 
Xj, since lagrangians are local; only these need be examined); 



if AS < (the action is reduced) retain the new value for Xj , and proceed 
to the next site; 



• if AS > generate a random number 77 unformly distributed between 
and 1; retain the new value for Xj if exp(— AS) > rj, otherwise restore 
the old value; proceed to the next site. 

An implementation of this algorithm, in the Python computer language 5 , is 
shown in Fig. 2. The code examples and Python are discussed in the Appendix. 
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def update (x) : 

for j in range(0,N): 

old_x = x[j] # save original value 

old_Sj = S(j ,x) 

x[j] = x[j] + uniform (-eps , eps) # update x[j] 

dS = S(j,x) - old_Sj # change in action 

if dS>0 and exp(-dS)<unif orm(0, 1) : 

x[j] = old_x # restore old value 

def S(j,x): # harm. osc. S 

jp = (j+l)°/„N # next site 

jm = (j-l)°/.N # previous site 

return a*x[j]**2/2 + x [j] * (x [j] -x [jp] -x [jm] )/a 

Figure 2: Python code for one Metropolis update of path {xj, j = . . . N — 1}. The path 
is stored in array x[j]. Function S(j,x) returns the value of the part of the action that 
depends on Xj. Function unif orm(a,b) returns a random number between a and b. A sample 
S(j ,x) is shown, for a harmonic oscillator with xn = xq. 



There are two important details concerning the tuning and use of this al- 
gorithm. First, in general some or many of the x/s will be the same in two 
successive random paths. The amount of such overlap is determined by the pa- 
rameter e: when e is very large, changes in the usually large and most 
will be rejected; when e is very small, changes are small and most are accepted, 
but the new Xj's will be almost equal to the old ones. Neither extreme is de- 
sirable since each leads to very small changes in x, thereby slowing down the 
numerical exploration of the space of all important paths. Typically e should 
be tuned so that 40%-60% of the Xj 's are changed on each pass (or "sweep" ) 
through the lattice. Then e is of order the typical quantum fluctuations ex- 
pected in the theory. Whatever the e, however, successive paths are going to 
be quite similar (that is "highly correlated" ) and so contain rather similar in- 
formation about the theory. Thus when we accumulate random paths x^ for 
our Monte Carlo estimates we should keep only every iV cor -th path; the inter- 
vening sweeps erase correlations, giving us configurations that are statistically 
independent. The optimal value for N COI depends upon the theory, and can 
be found by experimentation. It also depends on the lattice spacing a, going 
roughly as 

iV cor cx \. (35) 
Other algorithms exist for which iV cor grows only as 1/a when a is reduced, 
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but since our interest is in large a's we will not discuss these further. 

The second detail concerns the procedure for starting the algorithm. The 
very first configuration used to seed the whole process is usually fairly atyp- 
ical. Consequently we should discard some number of configurations at the 
beginning, before starting to collect x^'s. Discarding 5N COI to 10iV C or config- 
urations is usually adequate. This is called "thermalizing the lattice." 

To summarize, a computer code for a complete Monte Carlo calculation of 
((r[x])) for some function T[x] of a path x consists of the following steps: 

• initialize the path, for example, by setting all xj's to zero; 

• update the path 5-/Vc 0r -l(W C or times to thermalize it; 

• update the path 7V cor times, then compute T[x] and save it; repeat 
7V c f times. 

• average the N c { values of T[x] saved in the previous step to obtain a 
Monte Carlo estimate T for ({T[x})). 

A Python implementation of this procedure is shown in Fig. 3. 

Exercise: Write a computer program to implement the Metropolis Monte Carlo 
algorithm for a one dimensional harmonic oscillator (Eq. (18)), and calculate 

i 

for all t = 0, a, 2a . . . (N — l)a; that is calculate 

3 

for n = . . . N — 1. The (j + n)modiV in this last expression reflects the periodic 
boundary conditions. Try N — 20 lattice sites with lattice spacing a = 1/2, 
and set e = 1.4 and iV cor = 20. Try N c( 's of 25, 100, 1000 and 10000. Use the 
results to compute the excitation energy from 

AE n = log(G„/G„+i) n ^? c (Si - E Q )a (38) 

Try this for the harmonic oscillator potential and also for anharmonic poten- 
tials. Vary the various parameters. 

My results for the harmonic oscillator potential with N = 1000 configurations 
are shown in Fig. 4. These results required less than a minute of personal 
computer time. 
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def compute_G(x,n) : 
g = 

for j in range(0,N): 

g = g + x[j]*x[(j+n)°/.N] 
return g/N 

def MCaverage(x,G) : 

for j in range(O.N): # initialize x 

x[j] = 

for j in range (0, 5*N_cor) : # thermalize x 

update (x) 

for alpha in range (0 ,N_cf) : # loop on random paths 

for j in range (0 ,N_cor) : 

update (x) 
for n in range(0,N): 

G[alpha][n] = compute_G(x,n) 
for n in range(0,N): # compute MC averages 

avg_G = 

for alpha in range(0,N_cf ) : 

avg_G = avg_G + G [alpha] [n] 
avg_G = avg_G/N_cf 
print "G(7.d) = 7.g" 7„ (n,avg_G) 

Figure 3: Sample Python code for a Monte Carlo evaluation of of G(t) (Eq. (36)). Function 
computeG(x,t) computes G(t) for a given path x. Function MCaverage(x,G) computes the 
Monte Carlo average over random paths x. The results for path x' a ' are stored in the 
array G [alpha] [t] , and the averages are computed and printed. Function update (x) does 
one Metropolis sweep through the lattice (see Fig. 2). 
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AE(t) 




Figure 4: Monte Carlo values AE(t) = log(G(t)/G(t+a))/a plotted versus t for an harmonic 
oscillator. The exact asymptotic result, AE(oo) = 1, is indicated by a line. Results are for 
a one dimensional lattice with N = 20 sites, lattice spacing a = 1/2, and N c { = 1000 
configurations, keeping configurations only every N COT = 20 sweeps. The Metropolis step 
size e was 1.4, resulting in a Metropolis acceptance ratio of 0.5. 

Exercise: Redo the previous exercise but propagator 

Here we use x 3 rather than x to create and destroy the excited stated; that 
is, we use x 3 rather than x as the "source" and the "sink". Note that AE(t) 
converges to the same result, but only at much larger t's than before. Different 
sources and sinks often lead to different asymptotic behavior. Choices that 
result in fast convergence as t increases are usually preferable because statistical 
errors are smaller at smaller t's. Compare your best estimate of the asymptotic 
value obtained using x with that obtained using x. 

Note also that AE(t) approaches its its asymptotic value from above. Prove 
that this must be true in general, provided source and sink are the same oper- 
ator, using Eq. (23). This result is useful because it implies that each AE(t) 
gives a rigorous upper bound on the asymptotic value, even at small t's before 
convergence. 

My results with an x 3 source and sink are shown in Fig. 5. 
2.3 Statistical Errors 

A Monte Carlo estimate T of some expectation value ((r)) is never exact; there 
are always statistical errors that vanish only in the limit where infinitely many 
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Figure 5: Monte Carlo values AE(t) = log(G(t)/G(t+a))/a plotted versus t for an harmonic 
oscillator using x 3 as the source and sink. Parameters are the same as used to generate Fig. 4. 
The energies take longer to reach their asymptotic value. 

configurations are employed (N c { — > oo). An important part of any Monte 
Carlo analysis is the estimation of these statistical errors. There is a simple 
but very powerful method, called the "statistical bootstrap," for making such 
estimates. 

In the previous exercises, for example, we assemble an "ensemble" of mea- 
surements of the propagator G^ a \ one for each configuration x^ a \ These are 
averaged to obtain G, and, from it, an estimate for AE n (Eq. (38)). An obvious 
way to check the statistical errors on this estimate for AE n is to redo the whole 
calculation, say, 100 times, each time with different random numbers to gener- 
ate different random paths. With 100 copies of the entire calculation, we could 
analyze the distribution of the 100 random AE n 's obtained, and deduce the 
statistical uncertainty in our original estimate. This, however, is exceedingly 
expensive in computer time. The bootstrap procedure provides new, almost 
zero-cost random ensembles of measurements by synthesizing them from the 
original ensemble of N c f measurements. 

Given an ensemble {G^ a \ a = 1 . . . iV c f} of Monte Carlo measurements, we 
assemble a "bootstrap copy" of that ensemble by selecting G^'s at random 
from the original ensemble, taking N c f in all while allowing duplications and 
omissions. The resulting ensemble of G's might have two or three copies of 
some G^ 's, and no copies of others. This new ensemble can be averaged and a 
new estimate obtained for AE n . This procedure can be repeated to generated 
as many bootstrap copies of the original ensemble as we wish, and from each 
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def bootstrap (G) : 
N_cf = len(G) 

G_bootstrap = [] # new ensemble 

for i in range (0 ,N_cf) : 

alpha = int (unif orm(0 ,N_cf ) ) # choose random config 
G_bootstrap . append (G [alpha] ) # keep G [alpha] 

return G_bootstrap 

Figure 6: Sample Python code for producing a bootstrap copy of an ensemble of measure- 
ments G. The original ensemble consists of individual measurements G [alpha] , one for each 
configuration. The function bootstrap (G) returns a single bootstrap copy of ensemble G, con- 
sisting of N_cf measurements. Function unif orm(a,b) returns a random number between a 
and b. 



we can generate a new estimate for AE n . The distribution of these AE n 's 
approximates the distribution of AE n ''s that would have been obtained from 
the original Monte Carlo, and so can be used to estimate the statistical error 
in our original estimate. A Python implementation of this procedure is shown 
in Fig. 6. 

Another useful procedure related to statistical errors is "binning." At the 
end of a large simulation we might have 100's or even 100,000's of configu- 
rations x^ a \ and for each a set of measurements like G^ a \ our propagator. 
The measurements will inevitably be averaged, but we want to save the sepa- 
rate G^'s for making bootstrap error estimates and the like. We can save a 
lot of disk space, RAM, and CPU time by partially averaging or binning the 
measurements: For example, instead of storing each of 

G (l) G (2) G (3) G (4) G (5)_ (4Q) 

we might instead store 

,(D _ G« + G< 2 ) + G( 3 ) + GW 



4 

(2 ) _ G< 5 ) + G< 6 ) + G( ? ) + G® 



(41) 



The G^'s are far less numerous but have the same average, standard devia- 



tion, and other statistical properties as the original set. Typically the bin size 
is adjusted so that there are only 
this procedure is shown in Fig. 7. 



is adjusted so that there are only 50-100 G^'s. A Python implementation of 
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def bin(G,binsize) : 

G_binned = [] # binned ensemble 

for i in range (0,len(G) ,binsize) : # loop on bins 
G_avg = 

for j in range (0 ,binsize) : # loop on bin elements 

G_avg = G_avg + G[i+j] 
G_binned. append(G_avg/binsize) # keep bin avg 
return G_binned 

Figure 7: Sample Python code for producing a binned copy of an ensemble of measure- 
ments G. The original ensemble consists of individual measurements G [alpha] , one for each 
configuration. The function bin(G.binsize) bins the ensemble into bins of size binsize, 
averages the G's within each bin, and returns an ensemble consisting of the averages. 



Binning has an important side effect: it reduces or can even remove the 
effects of correlations between different configurations. If, when generating 
configurations, N COI is too small, successive Monte Carlo estimates are statis- 
tically correlated. This leads to error estimates, using Eq. (33), that can be 
much smaller than the true errors — a very bad situation. If, however, the 
Monte Carlo estimates are binned with sufficiently large bins, the majority of 
estimates in one bin will be uncorrelated from the majority in adjacent bins. 
Consequently the bin averages will be uncorrelated, and standard statistical 
formulas, like Eq. (33), are reliable. To determine the bin size required to 
remove correlations, first bin the measurements and compute the errors. Then 
rebin with double the bin size and recompute the errors. There are no corre- 
lations if the statistical errors are roughly independent of bin size; if, on the 
other hand, the statistical errors grow substantially with bin size (for example, 
proportional to the square root of the bin size), then there are strong cor- 
relations between bins. Continuing doubling the bin size until the statistical 
errors stop growing. Note that measurements of different physical quantities 
decorrclatc at different rates; different things may require different bin sizes. 

Exercise: Rerun your Metropolis simulation of the harmonic oscillator with N COI — 
1. Do several different runs and compare your results. Do they agree within 
statistical errors? Try binning the results from each of the runs in bins of 20 and 
recompute the statistical errors. Verify that different runs now agree within 
the errors computed from the binned results. 

In Fig. 8 I show results from N = 1000 configurations using N COI — 1. Compare 
this with results in Fig. 4, which come from the same number of configurations 
but with iV cor = 20. The errorbars in the N cor — 1 plot are obviously unreliable. 



14 



AE(t) 




Figure 8: Monte Carlo values AE(t) = log(G(t)/G(t + a))/a plotted versus t for an harmonic 
oscillator, as in Fig. 4 but with N CO r = 1. The errorbars are unreliable. 



3 Field Theory on a Lattice 

3.1 From Quantum Mechanics to Field Theory 

Field theories of the sort we are interested in have lagrangian formulations 
and so can be quantized immediately using path integrals. The procedure is 
precisely analogous to what we do in the previous section when quantizing 
the harmonic oscillator. The analogues of the coordinates x(t) in quantum 
mechanics are just the fields 4>(x) or A^(x) where x = (t,x) is a space-time 
point. Indeed our quantum mechanical examples can be thought of as field 
theory examples in spatial and 1 temporal dimension: x(t) — ► <f>{t) — > <p(x). 
The analogue of the ground state in quantum field theory is the vacuum state, 
|0), while the analogues of the excited states, created when <p(x) or cf> 3 or . . . acts 
on |0), correspond to states with one or more particles create in the vacuum. 
In the lattice approximation both space and time are discrete: 
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The nodes or "sites" of the grid are separated by lattice spacing a, and the 
length of a side of the grid is L; the lines joining adjacent sites are called 
"links." The quantum field is specified by its values at the grid sites: a con- 
figuration is describe by the set of numbers {<p(xj), Vxj-e grid}. The path 
integral generalizes in the obvious fashion: 

«rM»^ / c- s Mm Y[d^( Xj ) (42) 

Zjegrid 

where 

Z = J e~ SM Y[d<f>(xj). (43) 

Here the action S[<f>] is the continuum action with spatial and temporal deriva- 
tives replaced by differences between field values at the grid sites. We study 
excitations of the field theory using operators like 

r(t) = J= (44) 

Xj 

where the sum over the N spatial xj's enforces zero three-momentum. Since 
the excitations correspond to particle creation, their energies are the energies 
of particles: for example, 

«r(t)r(o)» ^ | <o| r(o) \4>-.p = o) | 2 e - m *«, (45) 

where \<p : p = 0) is a one ^-particle state with zero three momentum, and 
is the mass of the <j) particle. 

Exercise: Show that 

<0| T(0) |0 : p = 0) = J^ (46) 
where Z2 is the wavefunction renormalization parameter for the <j> field. 

3.2 Coarse Lattices 

We have, through the lattice approximation to the path integral, turned the 
problem of solving a nonperturbative relativistic quantum field theory, once 
again, into a problem of numerical integration. This is a major development for 
theories, like QCD, where perturbation theory doesn't suffice (at low energies). 

Early enthusiasm for such an approach to QCD, back when QCD was first 
invented, quickly gave way to the grim realization that very large computers 
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would be needed to numerically integrate the path integral. In recent years, 
however, two developments have made QCD simulations far more accessible. 
One is that small computers have become much faster; the other is that QCD 
simulations have become much faster — 10 3 to 10 6 times faster. These devel- 
opments imply that the simplest QCD simulations can be done using no more 
than a single personal computer or even a laptop. 

What has changed to make QCD simulations faster? The cost of a QCD 
calculation is given roughly by the formula 

costwf-) (47) 

where the first factor is the number of lattice sites, while the second and 
third factors are due to "critical-slowing-down" of the algorithms used for the 
simulation. From this formula, the single most important determinant of the 
cost is the lattice spacing: the cost is proportional to 1/a 6 . This means that 
one wants to keep a as large as possible. Until recently it was thought that 
a < 0.05-0.1 fm would be essential for accurate QCD simulations. As we shall 
see a « 0.3-0.4 fm works quite well. Given that the cost varies as 1/a 6 , the 
coarser lattices should be 10 3 to 10 6 times cheaper to simulate. 

The size of the lattice spacing is limited by discretization errors. The 
challenge is to make the lattice spacing as large as possible while keeping the 
discretization errors of order, say, a few percent or less. These errors have two 
sources: first, the lattice forces us to use approximate derivatives, and, second, 
it imposes an ultraviolet cutoff. We consider each in turn. 

In the lattice approximation, we only know the fields at the lattice sites. 
Thus all derivatives in field equations, the action, and the like must be con- 
verted to differences. For example, the second derivative of field </> at some 
point Xj on the lattice is given approximately by 

^l=A( 2 )0(x,) + O( a 2 ) (48) 



where 



A (2) = <l>(x + a)-2<j>{x)+<j>{x-a) ^ (49) 



We generally want more accurate approximations for work on coarse lattices: 
for example, the approximation 

^j^ 1 = Ai 2 ) - g (Ai 2 >) 2 fa) + 0(a 4 ) (50) 
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is accurate to within a few percent even when acting on structures in <j>(x) that 
are only four or five lattice spacings across. With such precision one might 
expect that lattice spacings as large as a quarter the diameter of hadron, or 
about 0.4 fm, would still be quite useful. Our theories are quantum theories, 
however, and therefore there is a second important consideration. 

The shortest wavelength oscillation that can be modeled on a lattice is one 
with wavelength A m ; n = 2a; for example, the function 4>{x) = +1, — 1, +1 . . . 
for x = 0, a, 2a . . . oscillates with this wavelength. Thus gluons and quarks 
with momenta p = 2ir/\ larger than ir/a are excluded from the lattice the- 
ory by the lattice; that is, the lattice functions as an ultraviolet cutoff. In 
simple classical field theories this is often irrelevant: short-wavelength ultravi- 
olet modes are either unexcited or decouple from the long-wavelength infrared 
modes of interest. However, in a noisy nonlinear theory like an interacting 
quantum field theory ultraviolet modes strongly affect infrared modes. Thus 
we cannot simply discard all particles with momenta larger than ir/a; we must 
somehow mimic their effects on infrared states. This is done by changing or 
"renormalizing" the parameters in our discretized theory and by adding new 
local interactions. 

The new interactions complicate the improved discretizations discussed 
above. For example, an interacting scalar theory on the lattice would have a 
discretized kinetic lagrangian 

E ¥' d > "> E 5 tA ^+ «W(A£°)V) (51) 

where parameter c has two parts: —1/12 from numerical analysis (Eq. (50)), 
and an additional renormalization due to the cutoff. Typically the renormaliza- 
tion is completely context dependent — for example, it is different for QED and 
QCD, or for particles of different spin, and so on. It cannot be looked up in a 
numerical analysis book; rather, it must be computed using quantum field the- 
ory. In QCD these renormalizations can be computed using (weak-coupling) 
perturbation theory, since the renormalizations are due to QCD physics at 
large momenta, p > it /a, where the theory is perturbative: 

c = --^ + CiO^Tr/a) + c 2 a 2 s (Tr/a) H . (52) 

This is true, that is, provided the lattice spacing is small enough that mo- 
mentum ir/a is perturbative. Work in continuum QCD suggests that lattice 
spacings of 0.5 fm or smaller should suffice, but, until recently, lattice simu- 
lations seemed to suggest that perturbation theory only started to work for 
lattice spacings smaller than 0.05-0.1 fm. 
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Exercise: Our action for one-dimensional quantum mechanics, Eq. (4), can be 
rewritten 

dt [-±mx{t)x(t) + V(x(t))] , (53) 

by integrating by parts (assuming x(ti) = x(t{) = x). To discretize we replace 

x( tj ) - A^x 3 ee ^-' 2 2 +X] -\ (54) 

where the Xj's are periodic (xo = xn, £-i = xn-i, and so on); this gives the 
same lattice action we used earlier. We can improve the discretization, however, 
by using the corrected approximation, Eq. (50), for the second derivative: 



S ilnp [x] = J2 a [-\ mx i ( AW - « 2 (A (2) ) 2 /12) Xj + aV{ Xj )] (55) 

3=0 

Modify your Monte Carlo code for the harmonic oscillator to include the cor- 
rection term. Run high-statistics comparisons (N c f — 10 4 or 10 5 ) with and 
without the correction term. 

Exercise: Discretized, euclidean classical equations of motion can be derived from 
the actions in the previous exercise by setting 



dS[x] 
dxi 



0. 



Using the improved action, for example, we obtain 

v (2) 2^(2)^2,^ _ _ dV( Xj 



m(A (2) -a 2 (A (2) ) 2 /12); 



dxj 



(56) 



(57) 



Setting 

V(x) = ^muigx 2 , (58) 

find solutions of the form Xj = exp(— uitj), the euclidean-time version of an 
oscillatory solution with frequency u>. Show that the frequency is given by 



2 2 



(59) 



for the unimproved action. The (acoo) 2 correction is the error caused by the 
finite lattice spacing. 

Repeat the exercise for the improved action and show that it has two solutions. 
One, 

(atuo) 4 



2 2 



1 + ■ 



90 



+ 0((au>f) 



(60) 
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Figure 9: Monte Carlo values AE(t) = log(G(t)/G(t+a))/a plotted versus t for an harmonic 
oscillator, as in Fig. 4 but with an improved action. The energies approach their asymptotic 
value from below. 



is an improved version of the previous result; its errors are fourth order in aula 
rather than second order. The other solution, however, is 



It corresponds to a new oscillation mode that does not appear in the continuum; 
it is an artifact of the improved lattice theory. This new mode is sometimes 
called a "numerical ghost." In a quantum field theory it would be a new, very 
massive particle (m oc 1/a). 

Our lattice theory was designed to be accurate for low-energies, and so we 
should not be surprised when unphysical modes appear at high energies. These 
ghost modes, being high-energy, typically decouple from low-energy physics 
and so can usually be ignored. However, they can have one unfortunate effect 
on the numerical analysis. Returning to the previous exercise, note that the 
AE n 's for the improved action are below the asymptotic result when n is small 
(see Fig. 9), in contradiction of the general result discussed in Section 2. The 
general result ignored the possibility that spurious states might be induced 
by the numerical analysis. Here these states have negative norm (impossible 
for real quantum states), which is why the energies rise from below. This 
is incovenient because it means that the AE n 's cannot be used to rigorously 
bound the true answer — they may be either above or below it — unlike the 
case for the unimproved action, where they must always be above. 
Ghost modes always arise when improved discretizations are used for temporal 
derivatives. There is a trick, however, for correcting temporal difference oper- 




(61) 
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ators that avoids extra states. This is to change integration variables in the 
path integral: for the harmonic oscillator we replace Xj — ► Xj where 

Xj — Xj + Sctj (62) 

and 

S£j = £1 a 2 A^Xj + £2 a 2 L0 2 Xj. (63) 
Substituting this into the action we obtain 

S[x] = S[x + Sx] 

ta dS[x] 



= S[x] (64) 

Find values for the such that the improved harmonic-oscillator action, in 
terms of x, is 

JV-l 

§im P [x] = ^2 a [-\ mS: 3^ {2) Xj + Vim P (a;j)] ( 65 ) 

3=0 

where 

V imp (x) = \mio 2 x 2 ^1 + { -^f-^j (66) 

and all 0(a 4 ) terms are ignored. This action has no a 2 errors, but also has 
no ghosts. The value of the path integral is not changed by a simple change 
of variables (provided that the jacobian is included — show that in this case it 
has no effect on expectation values). Rerun your numerical tests from the last 
exercise using this action. 

Exercise: The field transformation trick in the previous exercise is particularly sim- 
ple for the harmonic oscillator. Generalize the trick for the case of an anhar- 
monic oscillator with, for example, 

V(x) = imwoi 2 (l + cmLdox 2 ). (67) 

where c is a dimensionless parameter. Try a variable change with 

5xj = £i a 2 A^Xj + £2 a 2 u 2 Xj + £3 a 2 mu!ox^. (68) 



The resulting action is as above but with a new V\ 



mp ■ 



t> /~\ mui 2 _ 2 _2n a 2 mcjQ /_ _ 3 .,2 
Vimp(x) = — — x {1 + cmwox )H — — [x + 2cmaj x J 

3 

- aSv{£) + y Sv(£) 2 (69) 
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where the terms involving 

5v{x) = cmujox 2 /^ (70) 

are due to the jacobian (from the change of integration variables; the jacobian 
matters here because the change is nonlinear). Test this improved lattice action 
against the original unimproved action using m = u>o = 1 and c = 2. (The 
asymptotic value for AE is 1.933 with these parameters.) 
One might expect errors of order a 4 with the improved action. However renor- 
malization effects arise when interactions are anharmonic. In particular the 
coefficients of the x 2 and x 4 interactions are renormalized away from their 
naive values. Such effects enter at the same order as the corrections from the 
jacobian. In our one-dimensional theory, unlike in QCD, these corrections van- 
ish like powers of a when a — > 0. The leading such correction is an O(a) shift in 
the coefficient of the x 2 potential. Run high-precision simulations at a = 1/2 
and a = 1/4 to compute the O(a) error due to renormalization. Try adjusting 
the coefficient of x 2 to remove this error. (Alternatively one could try using 
perturbation theory to compute the shift needed to eliminate the error.) 

3.3 Perturbation Theory and Tadpole Improvement® 

Improved discretizations and large lattice spacings are old ideas, pioneered by 
Wilson, Symanzik and others . However, perturbation theory is essential; the 
lattice spacing a must be small enough so that pxur/a QCD is perturbative. 
This was the requirement that drove lattice QCD towards very costly simu- 
lations with tiny lattice spacings. Traditional perturbation theory for lattice 
QCD begins to fail at distances of order 1/20 to 1/10 fm, and therefore lattice 
spacings must be at least this small before improved actions are useful. This 
seems very odd since phenomenological applications of continuum perturbative 
QCD suggest that perturbation theory works well down to energies of order 
1 GeV, which corresponds to a lattice spacing of 0.6 fm. The breakthrough, 
in the early 1990's, was the discovery of a trivial modification of lattice QCD, 
called "tadpole improvement," that allows perturbation theory to work even 
at distances as large as 1/2 fm 6 ' 8 ' 9 . 

One can readily derive Feynman diagram rules for lattice QCD using the 
same techniques as in the continuum, but applied to the lattice lagrangian 11 . 
The particle propagators and interaction vertices are usually complicated func- 
tions of the momenta that become identical to their continuum analogues in 
the low-momentum limit. All loop momenta are cut off at p^ = ±w/a. 

Testing perturbation theory is also straightforward. One designs short- 
distance quantities that can be computed easily in a simulation (i.e., in a 

"This section is based upon work with Paul Mackenzie that is described in 6 . 
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Monte Carlo evaluation of the lattice path integral). The Monte Carlo gives 
the exact value which can then be compared with the pcrturbative expansion 
for the same quantity. An example of such a quantity is the expectation value 
of the Wilson loop operator, 

W{C) = (0|±Rc TrPe" is i A ' dx |0), (71) 

where A is the QCD vector potential, V denotes path ordering, and C is any 
small, closed path or loop on the lattice. W(C) is perturbative for sufficiently 
small loops C. We can test the utility of perturbation theory over any range 
of distances by varying the loop size while comparing numerical Monte Carlo 
results for W(C) with perturbation theory. 

Fig. 10 illustrates the highly unsatisfactory state of traditional lattice-QCD 
perturbation theory. It shows the "Creutz ratio" of 2a x 2a, 2a x a and a x a 
Wilson loops, 

_ f W(2ax2a)W(axa) \ 
X2 < 2 = ~H W>(2axa) )> (72) 

plotted versus the size 2a of the largest loop. Traditional perturbation theory 
(dotted lines) underestimates the exact result by factors of three or four for 
loops of order 1/2 fm; only when the loops are smaller than 1/20 fm does 
perturbation theory begin to give accurate results. 

The problem with traditional lattice-QCD perturbation theory is that the 
coupling it uses is much too small. The standard practice was to express 
perturbative expansions of short-distance lattice quantities in terms of the 
bare coupling a\ at used in the lattice lagrangian. This practice followed from 
the notion that the bare coupling in a cutoff theory is approximately equal to 
the running coupling evaluated at the cutoff scale, here a s (ir /a), and therefore 
that it is the appropriate coupling for quantities dominated by momenta near 
the cutoff. In fact the bare coupling in traditional lattice QCD is much smaller 
than true effective coupling at large lattice spacings: for example, 

«iat = a v (n/a) - 4.7 ay H (73) 

< \a v {^/a) fora>.lfm (74) 

where av{q) is a continuum coupling defined by the static-quark potential, 

^qq(^-4^ f ^M. (75) 

Consequently ai at expansions, though formally correct, badly underestimate 
perturbative effects, and converge poorly. 
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Figure 10: The X22 Creutz ratio of Wilson loops versus loop size. Results from Monte 
Carlo simulations (exact), and from tadpole-improved (new) and traditional (old) lattice 
perturbation theory are shown. 
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The anomalously small bare coupling in the traditional lattice theory is a 
symptom of the "tadpole problem" . As we discuss later, all gluonic operators 
in lattice QCD are built from the link operator 

[/■„ (*) ee 7>e -< XT +0 " « e -^ A " (76) 

rather than from the vector potential A^. Thus, for example, the leading term 
in the lagrangian that couples quarks and gluons is ipU^jftip/a. Such a term 
contains the usual tpgA ■ jip vertex, but, in addition, it contains vertices with 
any number of additional powers of agA^. These extra vertices are irrelevant 
for classical fields since they are suppressed by powers of the lattice spacing. 
For quantum fields, however, the situation is quite different since pairs of A^'s, 
if contracted with each other, generate ultraviolet divergent factors of 1/a 2 
that precisely cancel the extra o's. Consequently the contributions generated 
by the extra vertices are suppressed by powers of g 2 (not a), and turn out to 
be uncomfortably large. These are the tadpole contributions. 

The tadpoles result in large renormalizations — often as large as a factor of 
two or three — that spoil naive perturbation theory, and with it our intuition 
about the connection between lattice operators and the continuum. However 
tadpole contributions are gencrically process independent and so it is possible 
to measure their contribution in one quantity and then correct for them in all 
other quantities. 

The simplest way to do this is to cancel them out. The mean value uo 
of | Re Tr consists of only tadpoles and so we can largely cancel the tadpole 
contributions by dividing every link operator by uq. That is, in every lattice 
operator we replace 

U,(x) - (77 ) 

where u is computed numerically in a simulation. 

The uo's cancel tadpole contributions, making lattice operators and per- 
turbation theory far more continuum-like in their behavior. Thus, for example, 
the only change in the standard gluon action when it is tadpole-improved is 
that the new bare coupling an is enhanced by a factor of 1 /uq relative to the 
coupling ai a t in the unimproved theory: 

an = (78) 



Since Uq < .6 when a> .lfm, the tadpole-improved coupling is typically more 
than twice as large for coarse lattices. Expressing axi in terms of the continuum 
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coupling ay, we find that now our intuition is satisfied: 



«ti = a v (n/a) — .5 a v + • • • 
w ay(7r/a). 



(79) 
(80) 



Perturbation theory for the Creutz ratio Eq. (72) converges rapidly to the 
correct answer when it is reexpressed in terms of an- An even better result is 
obtained if the expansion is reexpressed as a series in a coupling constant de- 
fined in terms of a physical quantity like the static-quark potential, where that 
coupling constant is measured in a simulation. By measuring the coupling we 
automatically include any large renormalizations of the coupling due to tad- 
poles. It is important that the scale q* at which the running coupling constant 
is evaluated be chosen appropriately for the quantity being studied 6 ' 8 . When 
these refinements are added, perturbation theory is dramatically improved, 
and, as illustrated in Fig. 10, is still quite accurate for loops as large as 1/2 fm. 

This same conclusion follows from Fig. 11 which shows the value of the 
bare quark mass needed to obtain zero-mass pions using Wilson's lattice action 
for quarks. This quantity diverges linearly as the lattice spacing vanishes, and 
so should be quite perturbative. Here we see dramatic improvements as the 
tadpoles are removed first from the gluon action, through use of an improved 
coupling, and then also from the quark action. 

The Creutz ratio and the critical quark mass are both very similar to the 
couplings we need to compute for improved lagrangians. Tadpole improvement 
has been very successful in a wide range of applications. 

Asymptotic freedom implies that short-distance QCD is simple (pertur- 
bative) while long-distance QCD is difficult (nonperturbative). The lattice 
separates short from long distances, allowing us to exploit this dichotomy to 
create highly efficient algorithms for solving the entire theory: p>ir/a QCD is 
included via corrections SC to the lattice lagrangian that are computed using 
perturbation theory; p<ir/a QCD is handled nonperturbatively using Monte 
Carlo integration. Thus, while we wish to make the lattice spacing a as large 
as possible, we are constrained by two requirements. First a must be suffi- 
ciently small that our finite-difference approximations for derivatives in the 
lagrangian and field equations are sufficiently accurate. Second a must be suf- 
ficiently small that tt / a is a perturbative momentum. Numerical experiments 
indicate that both constraints can be satisfied when a w 1/2 fm or smaller, 
provided all lattice operators are tadpole improved. 
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Figure 11: The critical bare quark mass for Wilson's lattice quark action versus lattice 
spacing. Monte Carlo data points are compared with perturbation theories in a theory with 
no tadpole improvement (T.I), tadpole-improved gluon dynamics, and tadpole-improved 
quark and gluon dynamics. 
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4 QCD on a Lattice 

4-1 Classical Gluons 10 

The continuum action for QCD is 

S = f Jx^TrF^x) (81) 

where 

F„„ = d^A u - d v A„ + iglA^, A,} (82) 

is the field tensor, a traceless 3x3 hermitian matrix. The defining characteristic 
of the theory is its invariance with respect to gauge transformations where 

F„„ -» Q(x) Q(z)t (83) 

and Q(x) is an arbitrary x-dependent SU3 matrix. 

The standard discretization of this theory seems perverse at first sight. 
Rather than specifying the gauge field by the values of A fl (x) at the sites of 
the lattice, the field is specified by variables on the links joining the sites. In 
the classical theory, the "link variable" on the link joining a site at x to one 
at x + ajl is determined by the line integral of A^ along the link: 

/ r-x+afi \ 

U„(x) =Vexp -i / gA-dy) (84) 
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where the "P-operator path-orders the A^s along the integration path. We 
use f/^'s in place of A^'s on the lattice, because it is impossible to formulate a 
lattice version of QCD directly in terms of A^s that has exact gauge invariance. 
The U^'s, on the other hand, are SU3 matrices that transform very simply 
under a gauge transformation: 



t 



(85) 



This makes it easy to build a discrete theory with exact gauge invariance. 

A link variable U fJ- (x) is represented pictorially by a directed line from x 
to x + fi, where this line is the integration path for the line integral in the 
exponent of U^x): 



In the conjugate matrix U^{x) the direction of the line integral is nipped and 
so we represent U^(x) by a line going backwards from x + fi to x: 



A Wilson loop function, 



W(C) 



Ufa 



±TrPe~ 



i <j>^ gA-dx 



(86) 



for any closed path C built of links on the lattice can be computed from the 
path-ordered product of the U^'s and Ufa associated with each link. For 
example, if C is the loop 



then 



v 
t 



// 



W(C) = i Tr (U„(x)U„(x + o£) . . . Ufa)) 
28 



(87) 



Such quantities are obviously invariant under arbitrary gauge transforma- 
tions Eq. (85). 

One might wonder why we go to so much trouble to preserve gauge invari- 
ance when we quite willing give up Lorentz invariance, rotation invariance, etc. 
The reason is quite practical. With gauge invariance, the quark-gluon, threc- 
gluon, and four-gluon couplings in QCD are all equal, and the bare gluon mass 
is zero. Without gauge invariance, each of these couplings must be tuned inde- 
pendently and a gluon mass introduced if one is to recover QCD. Tuning this 
many parameters in a numerical simulation is very expensive. This is not much 
of a problem in the classical theory, where approximate gauge invariance keeps 
the couplings approximately equal; but it is serious in the quantum theory 
because quantum fluctuations (loop-effects) renormalize the various couplings 
differently in the absence of exact gauge invariance. So while it is quite possible 
to formulate lattice QCD directly in terms of A^'s, the resulting theory would 
have only approximate gauge invariance, and thus would be prohibitively ex- 
pensive to simulate. Symmetries like Lorentz invariance can be given up with 
little cost because the symmetries of the lattice, though far less restrictive, 
are still sufficient to prevent the introduction of new interactions with new 
couplings (at least to lowest order in a). 

We must now build a lattice lagrangian from the link operators. We require 
that the lagrangian be gauge invariant, local, and symmetric with respect to 
axis interchanges (which is all that is left of Lorentz invariance). The most 
local nontrivial gauge invariant object one can build from the link operators is 
the "plaqucttc operator," which involves the product of link variables around 
the smallest square at site x in the \iv plane: 

P^(x) = ± Re Tr (U^(x)U v (x + afj)U^(x + aft + av)Ut{x)) . (88) 

To see what this object is, consider evaluating the plaquette centered about a 
point xq for a very smooth weak classical field. In this limit, 

P»» « 1 (89) 

since 

Given that A^ is slowly varying, its value anywhere on the plaquette should 
be accurately specified by its value and derivatives at xq ■ Thus the corrections 
to Eq. (89) should be a polynomial in a with coefficients formed from gauge- 
invariant combinations of A f _ l (xo) and its derivatives: that is, 

P»v = 1 - ci a 4 Tr (gF^(xo)) 2 
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- c 2 a 6 Tr {gF^(x ){Vl + Vl)gF^{x )) 
+ 0(a 8 ) 



(91) 



where c\ and c 2 are constants, and D M is the gauge-covariant derivative. The 
leading correction is order o 4 because F 2 V is the lowest-dimension gauge- 
invariant combination of derivatives of A^, and it has dimension 4. (There 
are no F 3 terms because P^ u is invariant under — ► f/t or, equivalcntly, 
F -» -F.) 

It is straightforward to find the coefficients c\ and c 2 . We need only ex- 
amine terms in the expansion of that are quadratic in A^; the cubic and 
quartic parts of F 2 ^ then follow automatically, by gauge invariance. Because 
of the trace, the path ordering is irrelevant to this order. Thus 



/V = iReTrPe-^o^ 31 



= ± Re Tr 



1 



gA ■ dx 



gA-dx) +0(A 3 ) 



(92) 



where, by Stoke's Theorem 

r a/2 



cpA-dx = / dx^dx v [dpA^xa + x) - d u A^(x + . 

Jn J— a/2 



I 



/2 
a/2 

a/2 



dx^dx v [F M „(x ) + (x M D M + x^D^F^xq) 



« 2 ^(*o) + 24 + D 2 JF^(x ) + G(a 6 , A 2 ). 



(93) 



Thus a = 1/6 and c 2 = 1/72 in Eq. (91). 

The expansion in Eq. (91) is the classical analogue of an operator product 
expansion. Using this expansion, we find that the traditional "Wilson action" 
for gluons on a lattice, 



(94) 



where /3 — 6/g 2 , has the correct limit for small lattice spacing up to corrections 
of order a 2 : 



S wil = fd*x J2 h^F 2 ^ + ^TrF, v (Dl+Dl)F,„ + 



(95) 
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We can cancel the a 2 error in the Wilson action by adding other Wilson 
loops. For example, the 2a x a"rectangle operator" 



= | Re Tr 



(96) 



has expansion 



R^ = \-\o^ Tr (gF^) 2 - A a e Tr (gF^(4 D 2 + T^gF^) - • • • . (97) 

The mix of o 4 terms and a 6 terms in the rectangle is different from that in the 
plaquette. Therefore we can combine the two operators to obtain an improved 
classical lattice action that is accurate up to 0(a 4 ) 12 ' 13 : 



^classical 



X,fJ,>V 



5/v + R 



-I 



3 12 
d 4 x^iTr^ + 0(a 4 ). 



const 



(98) 
(99) 



This process is the analogue of improving the derivatives in discretizations of 
non-gauge theories. 6 

Exercise: Defining the "twisted-rectangle operator" 



T p „ = | Re Tr 



show that 



^classical 



-P { P »»+ i + 2 M }+const 



(100) 

(101) 
(102) 



This is an alternative to the improved gluon action derived in the previous 
exercise. 



6 An important step that I have not discussed is to show that the gluon action is positive for 
any configuration of link variables 13 . This guarantees that the classical ground state of the 
lattice action corresponds to F t _ lu =Q. 
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4-2 Quantum Gluons c 

In the previous section we derived improved classical actions for gluons that 
are accurate through order a 4 . We now turn these into quantum actions. The 
most important step is to tadpole improve the action by dividing each link 
operator by the mean link w : for example, the action built of plaquette 
and rectangle operators becomes 

The uo's cancel lattice tadpole contributions that otherwise would spoil weak- 
coupling perturbation theory in the lattice theory and undermine our procedure 
for improving the lattice discretization. Note that uo~3/4 when a = .4fm, and 
therefore the relative importance of the R^s is larger by a factor of l/u 2 ,~2 
than without tadpole improvement. Without tadpole improvement, we cancel 
only half of the a 2 error. 

The mean link u is computed numerically by guessing a value for use in 
the action, measuring the mean link in a simulation, and then readjusting the 
value used in the action accordingly. This tuning cycle converges rapidly to 
sclfconsistent values, and can be done very quickly using small lattice volumes. 
The uo's depend only on lattice spacing, and become equal to one as the lattice 
spacing vanishes. 

The expectation value of the link operator is gauge dependent. Thus to 
minimize gauge artifacts, uo is commonly defined as the Landau-gauge ex- 
pectation value, (0|i Tr {7 /x |0)lg- Landau gauge is the axis-symmetric gauge 
that maximizes uq, thereby minimizing the tadpole contribution; any tadpole 
contribution that is left in Landau gauge cannot be a gauge artifact. An alter- 
native procedure is to define u as the fourth root of the plaquette expectation 
value, 

- (OI/VIO} 1 / 4 . (104) 

This definition gives almost identical results and is more convenient for numer- 
ical work since gauge fixing is unnecessary. 

Tadpole improvement is the first step in a systematic procedure for improv- 
ing the action. The next step is to add in renormalizations due to contributions 
from k>w/a physics not already included in the tadpole improvement. These 
renormalizations induce a 2 a s (ir / a) corrections, 

6£ = a s r ia 2 Y, Tr(^D^) 

c This section is based on work with M. Alford, W. Dimm, G. Hockncy and P. Mackenzie that 
is described in 14 . 
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+ a s r 2 a 

Fya ) 

+ a s r 3 a 2 £ Trp^D^) 

+ (105) 

that must be removed. The last term is harmless; its coefficient can be set to 
zero by a change of field variable (in the path integral) of the form 

i4„-> A„ + a 2 a s f(a s ) ^D^. (106) 

Since changing integration variables does not change the value of an integral, 
such field transformations must leave the physics unchanged. d Operators that 
can be removed by a field transformation are called "redundant." The other 
corrections are removed by renormalizing the coefficient of the rectangle oper- 
ator R^ v in the action, and by adding an additional operator. One choice for 
the extra operator is 



C^ua = \ Re Tr 



(107) 



Then the action, correct up to 0(a 2 a 2 ,a 4 ), is 



15 



^-'Efi^-'^}^' £ < 108) 

where 

r g = 1 + .48a s (7r/a) (109) 

c g = .055a s (7r/a). (110) 

The coefficients r g and c g are computed by "matching" physical quantities, 
like low-energy scattering amplitudes, computed using perturbation theory in 
the lattice theory with the analogous quantity in the continuum theory. The 
lattice result depends upon r g and c g ; these parameters are tuned until the 
lattice amplitude agrees with the continuum amplitude to the order in a and a s 
required: 

7iat(r g ,C g ) = Tcontin- (HI) 

d One must, of course, include the jacobian for the transformation in the transformed path 
integral. This contributes only in one-loop order and higher; it has no effect on tree-level 
calculations. 
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Note that tadpole improvement has a big effect on these coefficients. Without 
tadpole improvement, r g = 1 + 2a s ; that is, the coefficient of the radiative 
correction is four times larger. Tadpole improvement automatically supplies 
75% of the one-loop contribution needed without improvement. Since a s K0.3, 
the unimproved expansion for r g is not particularly convergent. However, with 
tadpole improvement, the one-loop correction is only about 10-20% of r g . 
Indeed, for most current applications, one-loop corrections to tadpole-improved 
actions are negligible. 

Exercise: Show that the gauge that maximizes (0| | Tr £/ M |0) becomes Landau gauge 
(3 ■ ,4 = 0) in the a->0 limit. 

4-3 Monte Carlo Evaluation of Gluonic Path Integrals 

A computer code for the Monte Carlo evaluation of gluonic path integrals 
can be designed in close analogy with our code for one-dimensional quantum 
mechanics. The Metropolis algorithm for generating random configurations 
must be adapted to work with SU3 matrices. In our quantum mechanics ex- 
ample, the coordinate was updated by adding a random number. In QCD 
the gluon field is specified by link variables U^(x), which are exponentials of 
the fundamental field. Thus to update a link variable we must multiply by 
the exponential of a random field; that is we must multiply by a random SU3 
matrix M: 

U^MU^ (112) 

Typically the matrix M is chosen randomly from a set of 50 or 100 random 
SU3 matrices that is generated once, at the start of the simulation; the only 
restrictions on this set are that it be large enough so that products of the 
various M's cover the entire space of SU3 matrices, and that the inverse, M\ 
for each matrix M in the set also be included in the set. The M's can be 
generated by first creating a set of hermitian matrices H whose matrix elements 
are random numbers between —1 and 1. These are converted to SU3 matrices 
by forming 1 + ieH and unitarizing it. e Parameter e determines the size of the 
update; as before, it is adjusted so that roughly half of all trial updates are 
accepted. 

A second modification of the Metropolis algorithm that is useful for QCD is 
to update each link variable several times (rather than just once) before moving 
on to the next variable in a single sweep through the lattice. This allows the 

e To convert an arbitrary matrix M = (mi m,2 m^) into an SU3 matrix: first normalize the 
first column mi to unity, then make the second column m,2 orthogonal to the (new) first 
column and normalize it, and replace the third column by the cross product of the (new) 
first two columns. 
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link variable to come into statistical equilibrium with its immediate neighbors 
on the lattice. The additional cost for these extra updates is relatively small 
because standard gluon actions are linear in the link variables. Thus the part 
of the action that must be computed when updating a particular U^(x) can 
be written 

AS{x,n) = Re Tr {U^T^x)) (113) 

where T fJ _(x), which is a sum of products of the link variables, is independent 
of Uft(x). Computing T^(x) is the most expensive part of the Metropolis 
update, but it need be computed only once for each set of successive updates 
of Uft(x). Typically one does about 10 "hits" of the Metropolis algorithm 
before moving on the the next link variable. 

Exercise: Design a computer code for evaluating gluonic path integrals using the 
Metropolis algorithm. Do this first for the simplest lattice action, the Wilson 
action (Eq. (94)): 

sw^-^E^# ( 114 ) 

X,[l>h< 

Run simulations at (3 = /3/uq = 5.5, which corresponds to a lattice spacing of 
around 0.25 fm. The lattice volume should be of order 2 fm on a side for typical 
QCD simulations; use L/a = 8 points on a side in your simulation. Set the 
Metropolis step size e = 0.24 and omit 7V cor = 50 sweeps between Monte Carlo 
measurements. Compute Monte Carlo averages of a x a and a x 2a Wilson 
loops; you should obtain about 0.50 and 0.26, respectively. 
Also try the improved action, Eq. (103). Use f3 = 1.719 and uo = 0.797 to 
again obtain a « 0.25 fm 16 . The a x a and a x 2a Wilson loops have values 
of 0.54 and 0.28, respectively. (Wilson loops are unrenormalized and so these 
values need not agree with those from the Wilson action.) 

4-4 A First Simulation 

Perhaps the simplest physics calculation that involves just gluons is to compute 
the potential energy between a static quark and a static antiquark separated 
by a distance r. This "static potential" should be approximately Coulombic 
at short distances, but grow linearly at large distances, demonstrating quark 
confinement. It can be used in a Schrodinger equation to predict energy levels 
for the ip/J and T families of mesons. 

The euclidean Green's function or propagator G for a heavy nonrelativistic 
quark in a background gauge field satisfies the equation 

(d*-7^)g(x)=«5 4 (x) (115) 
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where D p = — igA^(x) is the gauge-covariant derivative. This equation is 
easily solved in the static-quark limit, where M — > oo, to obtain 



Goo(x,i) 



r Pd' i Jo 9 A o(x-,t)dt 



t 



<5 3 (x), 



which on the lattice becomes 

G oc (x, t) = E/ftx, t-a) E/ftx, t-2a) . . . C//(x, 0). 



(116) 



(117) 



Propagation of a static antiquark is described by G^. Therefore we obtain 
the static potential V(r), which is the energy of a static quark and antiquark 
a distance r apart, from expectation values ofrxt Wilson loops: 



W{r,t) = (0\±Tr 



|0) r 



(118) 



where for large t 



W(r,t) const e - y(r)t . 



(119) 



Thus, to calculate the static potential V(r), we compute W(r, t) for a variety 
of t's and then take the large-t limit 



W(r, t)/W(r,t + a) -» aV(r). 



(120) 



There is one modification of this procedure that greatly improves the 
results. This is to replace the spatial link matrices, in the r direction, by 
"smeared" link variables U^x): 



where 



%(x) = (l + ea 2 A( 2 ))" U^x) 



(121) 
(122) 



(2) 

and Ap is a discrctizcd, gauge-covariant derivative: 



^ ] U,{x) ee 



,2 „2 



(U p (x) U„(x + ap) UUx + ap)-2 u 2 U^x) 



+ U^(x - ap) Ufj,{x - ap) U p (x - ap + ap)) . (123) 
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The low-momentum components of the smeared link variable are unchanged 
by the smearing, up to corrections of order a 2 p 2 ; but at high momentum the 
smearing acts as an ultraviolet cutoff: 



(f + ea 2 A( 2 ))"^ (1 




(124) 



By suppressing high-momentum gluons in the initial and final states of our 
matrix element, we suppress contributions from gluonic excitations and thereby 
hasten the convergence of W(r, t)/W(r, t + a) to its asymptotic value. This 
significantly reduces the Monte Carlo statistical errors since these are generally 
much smaller at smaller t's. 

Exercise: Add measurement code for W(r, t) to your gluon Monte Carlo program. 
Run this code for the parameter sets given above and extract the static poten- 
tial. Try this with unsmeared spatial (7 M 's on the ends of the Wilson loop, and 
again with smeared links (try n = 4 smearings with e — 1/12). In the smeared 
case asymptotic results appear at the first or second time step 
My results for the Wilson and improved actions are shown in Fig. 12. Each run 
took about 3 hours on a 300 MHz personal computer (with 64MB of RAM). 
The potential was computed for quarks separated along the lattice axes, as 
well as for separations along diagonal directions on the lattice (using nonplanar 
Wilson loops). The plots clearly show the linear rise of the potential at large r, 
which results in quark confinement; and careful examination shows the onset 
of Coulombic behavior at small r, which is due to asymptotic freedom. The 
Monte Carlo results are compared in the figure with global fits to the form 



where a is the string tension. This parameterization works well for the range of 
r's we are considering (0.25 fm-1 fm). Note that the Monte Carlo data is not as 
smooth in the Wilson case as it is for the improved action. This is caused the 
0(a 2 ) errors in the Wilson action, coming from the second term in Eq. (95). 
This term breaks rotational invariance; it is the leading manifestation in the 
gluon dynamics of the cubic structure of our lattice. While the true potential is 
a function only of the magnitude of the quark separation, we expect differences 
on the lattice between potentials for separations that are along the grid axes 
and those for separations along diagonals. This is particularly apparent at 
r = 3a where there are two Monte Carlo points, one for separation (3a, 0, 0) 
and the other for (2a, 2a, a). The two points are clearly separated on the 
plot for the Wilson action. The a 2 errors are removed from the improved 
action; the two r = 3a points merge on the plot for that action. The difference 
a(V(2a, 2a, a)-V(3a, 0, 0)) is reduced from 0.065±0.007 with the Wilson action 
to 0.003 ± 0.006 with the improved action in these particular simulations. The 
a 2 errors are larger at smaller r's, but more difficult to quantify there. 



V(r) — or — b/r + c 



(125) 
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Figure 12: Static-quark potential computed using the Wilson gluon action and the 0(a 2 ) 
improved action. The dotted line in each case is a fit of the Monte Carlo results to crr—b/r+c. 
The lattice spacing a 0.25 fm in each case. Each plot required about 3 hours of computer 
time using a 300 MHz personal computer. The improved action gives a smoother curve. 
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5 Conclusions 

The QCD simulation discussed in the previous section is, of course, just a 
beginning. A reasonable next step would be to add heavy quarks to the simu- 
lation, using nonrelativistic QCD, and to compute spectra and wavefunctions 
for tjj/J and T mesons 10, 17, 18 . Far more costly to simulate, although no more 
subtle theoretically, are light quarks 10,16 . These are needed to analyze protons, 
neutrons and other standard hadrons. Finally, and most costly, one might 
include light-quark vacuum polarization, which is an essential step for high- 
precision calculations (better than 10-15%). 
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Appendix 

The sample simulation code in Section 2 is written in Python, a simple but 
powerful computer language that is freely available for just about any type of 
computer from the web site www.python.org. These simulations are most effi- 
cient when the Numeric package is used with Python; this is included with some 
Python distributions but must downloaded separately (from www.python.org) 
for others. To use the code given in the text, collect it together in a single file 
as follows: 

import Numeric 

from whrandom import uniform 
from math import * 

# ... code from text goes here 

# set parameters: 
N = 20 

N_cor = 20 
N_cf = 100 
a = 0.5 
eps = 1.4 

# create arrays: 

x = Numeric . zeros ( (N ,) .Numeric . Float) 

G = Numeric . zeros ( (N_cf ,N) .Numeric . Float) 
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# do the simulation: 
MCaverage(x,G) 



If the file is called simulation, py, it is run with the command python simulation 
To test the binning and bootstrap codes add the following to the the file: 



g = Numeric . asarray(G) 

return Numeric . absolute (avg(g**2) -avg(g) **2) **0 . 5 

print 'avg G\n',avg(G) 

print 'avg G (binned) \n' ,avg(bin(G, 4)) 

print 'avg G (bootstrap) \n' , avg (bootstrap (G) ) 

The average of the binned copy of G should be the same as the average of G 
itself; the average of the bootstrap copy should be different by an amount of 
order the Monte Carlo error. Compute averages for several bootstrap copies 
to get a good feel for the errors. 

Finally one wants to extract energies. This is done by adding code to 
compute AE(t): 

def deltaE(G): # Delta E(t) 

avgG = avg(G) 

adE = Numeric . log(Numeric. absolute (avgG [ : -1] /avgG[l :] )) 
return adE/a 

print 'Delta E\n' ,deltaE(G) 

print 'Delta E (bootstrap) \n' , deltaE(bootstrap(G) ) 

Again repeating the evaluation for 50 or 100 bootstrap copies of G gives an 
estimate of the statistical errors in the energies. Additional code can be added 
to evaluate standard deviations from these copies: 

def bootstrap_deltaE(G,nbstrap=100) : # Delta E + errors 



def avg(G) : 

return Numeric . sum (G) /len(G) 



# MC avg of G 



def sdev(G) : 



# std dev of G 



avgE = deltaE(G) 
bsE = [] 



# avg deltaE 
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for i in range (nbstrap) : # bs copies of deltaE 

g = bootstrap (G) 

bsE . append (deltaE (g) ) 
bsE = Numeric . array (bsE) 

sdevE = sdev(bsE) # spread of deltaE' s 

print "\n°/„2s °/ 10s °/„10s" % ("t" , "Delta E(t) ", "error") 
print 26*"-" 

for i in range(len(avgE)/2) : 

print "7„2d °/„10g °/„10g" 7, (i , avgE [i] , sdevE [i] ) 

bootstrap_deltaE(G) 

The entire program should take only a few seconds to run on a modern personal 
computer. 

This example almost completely ignores the powerful object-oriented fea- 
tures that distinguish Python from most other scripting languages. 
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